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Abstract 

We present a theory for the damping of layer-by-layer growth oscillations in molecular 
beam epitaxy. The surface becomes rough on distances larger than a layer coherence length 
which is substantially larger than the diffusion length. The damping time can be calcu- 
lated by a comparison of the competing roughening and smoothening mechanisms. The 
dependence on the growth conditions, temperature and deposition rate, is characterized by 
a power law. The theoretical results are confirmed by computer simulations. 
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I. INTRODUCTION 



Basic questions of crystal growth have been a source of inspiration in physics for many 
years 0. Several concepts of general theoretical importance were originally developed in 
this area, but a good deal of the fascination is also due to the technological relevance. 

Here we want to address layer-by-layer growth or Frank-van der Merwe growth, as it is 
also called . This is a growth mode in molecular beam epitaxy which allows well controlled 
manipulation of e.g. chemical composition in layers down to atomic thickness. It is therefore 
particularly suited for the fabrication of electronic devices. 

In the most elementary model of layer-by-layer growth atoms are deposited under high 
vacuum conditions onto a (100) or (111) surface, for instance. The adatoms diffuse on 
the surface until they meet one another to form dimers which then grow into islands of 
monoatomic height. Their edges capture most of the adatoms during the deposition of one 
monolayer. When the available island edges become less due to coalescence, the formation of 
dimers and islands in the next layer begins. Hence the density of atomic steps like all other 
quantities which are sensitive to the surface morphology oscillate. These oscillations are 
the real hallmark of layer-by-layer growth. They can be used to monitor the film thickness 
easily. 

In step flow growth of a vicinal surface, where the nucleation of islands on the terraces 
can be neglected, as well as for rough surfaces the oscillations are absent. In practice, one 
cannot prepare substrates without a miscut. It may happen to be very weak, but it cannot 
be controlled. Layer-by-layer growth with its oscillations can only be seen if some terraces 
between steps are wide enough to allow for the formation of islands. 

Generically, the oscillations are damped: Layer- by-layer growth is only a transient. The 
reason can be either the emergence of a smooth surface with step flow or of roughness [[|, 
two entirely different diagnoses for the same symptom. However, increasing the temperature 
the damping should become stronger in the first and weaker in the second case, allowing to 
discriminate the two. Roughness itself can have two different reasons. If interlayer transport 
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is inhibited by step edge barriers, adatoms accumulate on top of islands, and one obtains 
the growth instability predicted by J. Villain 0, demonstrated and analysed in computer 
simulations || and observed in many experiments 0. If no instability occurs, the surface 
may still roughen due to fluctuations in the beam intensity (shot noise). This is the case 
which will be discussed below in more detail. 

Layer-by- layer growth is also possible if one starts from a patterned substrate rather 
than a perfectly smooth one, and one may wish to optimize the growth conditions such that 
the pattern stays intact for several layers. This optimization process is nontrivial: On one 
hand, the islands nucleate with a characteristic distance I, which sets the spatial resolution 
for the pattern in subsequent layers. On the other hand, deposition, diffusion and island 
nucleation are stochastic processes and lead to the accumulation of errors in the pattern, 
the thicker the film grows. Often, reducing these errors implies an increase of I and hence a 
worse resolution of the pattern. 

In order to understand layer-by- layer growth, one therefore has to investigate spatial 
correlations, temporal correlations and finally the emergence of surface roughness, which 
determines, for how many layers the oscillations can be seen. 

The spatial correlations in layer-by- layer growth manifest themselves through the char- 
acteristic distance / of nucleation events within a layer mentioned above. The dependence 
of I on growth conditions has already been well studied |7Hl3|, and the key result needed 
later is 

I ex (D/F)\ (1) 

where D and F are the surface diffusion constant and the deposition rate, respectively. The 
exponent 7 depends on the adatom diffusion process and the (fractal) dimension of the 
islands. It also depends on whether or not desorption of adatoms or diffusion of dimers or 
larger clusters are negligible, as we assume in the simulations presented below. Finally, 7 
is a function of the critical island size i*, which is defined by the size i* + 1 of the smallest 
island, which is stable enough that it never decays before capturing the next adatom. 



The temporal correlations are interesting, because they determine how well a nucleation 
pattern is reproduced in subsequent layers. For example there is a high probability that 
the first island nucleates on top of the first one of the previous layer, because this is the 
area where adatoms accumulate most likely. In fact, it was shown that the autocorrelation 
function of the first nucleation events decays rather slowly with time like t -1 / 2 Of 



course, this power law is only valid as long as layer-by-layer growth persists. It will be cut 
off by the time t, which characterizes the damping of the oscillations. 

The main subject of this paper is the emergence of surface roughness responsible for the 
damping of the oscillations. The simulation results show that layer-by-layer growth goes on 
forever if the linear size of the system is smaller than a layer coherence length /. Up to this 
length the layers grow coherently, for larger distances they get out of phase. Remarkably, 
/ is much larger than the characteristic distance between islands, I. The surface becomes 
rough on scales larger than / rather than the diffusion length I. 

In order to study kinetic roughening one may average the film thickness over the distance 
I. Then one cannot resolve individual islands any more, but still sees the dephasing between 
layers. Phenomena on this scale can be described by continuum equations, which provide the 
most transparent theoretical framework in which to discuss the smoothening mechanisms 
competing with the shot noise. The layer coherence length I as well as the damping time 
t play an important role for kinetic roughening as natural cutoffs of the continuum growth 
equation at small length and time scales. This idea will be worked out in section ||. 



Recently, it has been discovered [15,12| that t depends on the growth conditions like 



Ft oc (D/F) 5 . (2) 

This law will be explained in section [TTI| based on the results of section 0. Both sections, 
1 1] and |TTT| rely on dimensional arguments. In section [TV] the key results of the preceding 
section are rederived avoiding dimensional arguments. Section [V] contains a discussion of 
crossovers due to different smoothening mechanisms, and section [VT] reviews the numerical 
results confirming the theoretical picture. 



II. A THEORY FOR THE LAYER COHERENCE LENGTH AND THE 

DAMPING TIME 



The transition from layer-by-layer growth with its oscillations to kinetic roughening 
happens at time t, when the film thickness varies over the distance / by about one atomic 
layer. For t > t one expects that the surface shows self affine scaling ]IE| : 



w{t) « a±(£(t)/fy with f (t) « I (t/t) 1/z . (3) 

Here, w is the root mean square variation of the film thickness, a± the thickness of one atomic 
layer, and £ the correlation length up to which the surface roughness has fully developed 
until time t. ( is the roughness exponent and z the dynamical exponent. 

t is the time at which a continuum description of kinetic roughening becomes appropriate. 
Whenever desorption and the formation of defects in the growing film can be neglected the 
equation of motion must have the form 

dth = -V-i+ri, (4) 

where h is the deviation of the film thickness from its average value and rj(x,t) denotes the 
shot noise with correlator 

d is the surface dimension. In the conserved KPZ (cKPZ) equation, which has been widely 
discussed (ffl|T7|-[I!| in the context of molecular beam epitaxy, the adatom current has two 
terms, one driven by differences in the surface curvature and the second one by differences 
in the squared surface tilt: 

j = V[i^V 2 HA(V/i) 2 ]. (6) 

Eq. (|J) shows that the only characteristic length, time and height entering the description 
of the rough surface (coarse grained on scale /) are I, t and a±, respectively. Therefore K, 
A and JF must be functions of these three quantities. For example, A has the dimension 



length 4 height 1 time 1 . This implies that it must be the product of a dimensionless factor 
and l 4 /a±i. Similarly one obtains 

K ps o_lA ps Z~ 4 /t (7) 

According to (fD, 77 has the dimension height/time. Taking the dimensions of the (^-functions 
in (|^) (length - ^ and time -1 , respectively) into account, one finds that T is 

T ps a 2 J d /t, (8) 

up to a dimensionless factor. 



III. CONNECTION WITH SUBMONOLAYER PHYSICS 

In order to derive (@) from (|7|), flD one has to know, how K (or aj_A) and T depend on 
D and F. This question will be answered in the following. 

The physics of kinetic roughening should be determined by the same microscopic pro- 
cesses that are also responsible for the phenomena in the submonolayer regime. There, the 
characteristic time is the layer completion time, 

r = {Fa d )-\ (9) 

and there are two characteristic lengths, the diffusion length, I, and the lattice constant 
along the surface, a. Therefore, it must be possible to express K, X and T in terms of I, a, 
t and a±. 

The coefficients A and K characterize the morphology dependence of the nonequilibrium 
adatom density, which drives the surface current (see Section |TV| ). The most important 
morphological feature is the typical distance between islands. Therefore it is natural to 
assume that K and aj_X are only functions of I and r. The only dimensionally correct 



expressions are then EO] 



K ps a ± \ ps / 4 /r. (10) 



By contrast, the shot noise cannot depend on surface diffusion. Therefore, the diffusion 
length I cannot enter, and T /a\ can only be a function of a and r. The only dimensionally 
correct expression is then 



T « a\a d /T = F(a ± a d ) 2 (11) 



Comparing (0,0) with (|7|||) one finds that 

! 4 /t^l 4 /r and ! d /t^a d /r. (12) 

This, finally, leads to the central result of this paper, 

l/a « (l/af^ and Fa d t = t/r « {l/af d '^- d \ (13) 

Note in particular that the layer coherence length I I. With (|]) the exponent S defined 
in © is 

provided the cKPZ equation is the appropriate continuum equation for the growth process. 
At the upper critical dimensionality d — d c — 4 the scales / and t depend exponentially on 
I, while for d > d c the oscillations persist forever and the surface remains smooth. 

IV. THE ADATOM CURRENT AND THE SHOT NOISE RECONSIDERED 

In this section ([H]) and (pT|) will be rederived without using dimensional arguments. 
First we give a microscopic derivation of the nonlinear contribution to the adatom current 
(|) (see also |H|). 

It was proposed by Villain [|]] that in growth processes far from equilibrium, where local 
chemical potentials along the surface are ill defined, diffusion currents should be driven by 
gradients in the growth-induced, nonequilibrium adatom density n, 
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j = -Da L a d Vn. (15) 

The factor a d a± is the atomic volume and enters because @ expresses volume rather than 
mass conservation. 

On a singular surface the balance between deposition and capture of adatoms at steps 
leads to a stationary adatom density n = no of the order |K| 



n » (F/D)l 2 . (16) 

On a vicinal surface the adatom density is reduced due to the presence of additional steps; 
however this effect is felt only if the miscut m = |V/i| exceeds a±/l, in which case (16) is 
replaced by n « (F/D)(a±/m) 2 . In terms of a coarse grained description of the surface this 
implies that the local adatom density depends on the local miscut or surface tilt. A useful 



interpolation formula which connects the regimes m <C a±/l and m > ai/I is [20 



#) = l + {l v h/a±) 2 w ~ iF/D)l\Vh/a L f + ... (17) 



Inserting the leading quadratic term of this gradient expansion into ([151) , which is appropriate 
for describing long wavelength fluctuations around the singular orientation, we obtain 

j = VA(V/i) 2 (18) 

with A = Fa d l 4 /a±, which agrees with the result ( |I0"D of the previous section. 
Now we rederive ( |TT| ) for the shot noise, which can be written as 

l d /t = Fa 2d (19) 

by inserting (§) and (|j). This follows from considering the number of particles deposited 
during time i into area l d , Ft l d ± (Ft l d ) 1 ^ 2 . Each particle contributes a volume a d a±, so 
that the fluctuation of the film thickness over the distance I is 



w(t) w \IFt l d a d a ± /l d . (20) 
At t this should be the thickness of about one atomic layer, w(t) ~ aj_, which results in (|T9|). 



V. COMPETING MECHANISMS 



Whereas for our computer simulation model the theoretical arguments given in the pre- 
ceding sections are perfectly appropriate, the question arises how relevant these results are 
in general experimental situations. It has been argued that generically one should expect 
nonequilibrium contributions to the surface current which are driven by a height difference 
P JT7|j2^| . To leading order in a gradient expansion one gets an adatom current of the form 

j = -uVh. (21) 



Eq. (f|) with such a current is known as the Edwards- Wilkinson (EW) equation [^3[ . 

Tilt induced nonequilibrium surface currents originate from step edge barriers of Ehrlich- 
Schwoebel-type P4| , as well as kick-out or diffusion exchange processes at step edges. 



Whereas the latter two lead to a downhill current stabilizing the surface [y > 0), the 
former generates an uphill current [y < 0) and consequently an instability which will not be 
considered in the following. 

In the case of kick-out processes the coefficient v cannot depend on the diffusion length, 
because they are caused by deposition events in the immediate vicinity of a downward step. 
The only dimensionally correct expression is therefore 

u = a 2 /r = Fa d+2 . (22) 

The corresponding current is proportional to the local step density V/i and the deposition 
rate F. Such contributions will be absent in the simulations discussed in the next section. 
In general, the adatom current will contain the terms @ as well as (|2l|). The latter 



one dominates the surface roughness on large scales. Whether or not it influences the 
damping of the growth oscillations, however, depends on the crossover time t\ v from cKPZ- 
(A-dominated) at early to EW- (^-dominated) behaviour at late times. If the oscillations 
are damped out before the crossover takes place, the A-term determines the damping, hence 
the above result applies. Let t\ and i u denote the damping times if only the A- or the z/-term 
were present in the continuum equation of motion. Then ( |l3l) and ([14]) hold if 
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tx < t X u- (23) 

If, however, this is not the case, then t\ is replaced by £„, as long as no further terms in the 
continuum description provide further time scales. 

The crossover time t v \ is estimated in the following way: First we calculate the typical 
height fluctuation h u (t) after time t, if the A-term would be absent. Similarly, h\(t) is the 
fluctuation amplitude, if v = 0. Equating h x and h v then gives t\ v . By dimensional analysis 
one gets (see appendix) 

^ \ 4/(d+2) ^j-.^ (d+8)/(d+2) d ( l\ 16/(d+2) 



where (0-0) and (^) have been used to replace the parameters A, T and v. 
This has to be compared with (|T3"D, 

/ , v 4d/(4-d) 

« aia d ( ij . (25) 

For d < 2 the damping time is smaller or of equal order of magnitude as the crossover time 
t\ u . This implies that kick-out processes at step edges, although leading to an EW-term 
in the growth equation and hence modifying the later roughness, do not change our results 



(|T3|) for the layer coherence length and the damping time. 

However, if for example the sticking probability at an up step would be much smaller 
than at a down step (e.g. due to a step decoration by surfactant atoms p5|), one would 
expect a downhill current depending on I rather than the lattice constant a, i.e. with 

v w I 2 /t (26) 

instead of (|22|). In this case ( p4"| ) is replaced by 

/ j v -2d/(2+d) 

Tt Xv w a ± a d - , (27) 



a 



which is never larger than Tt\. The damping time should then be given by ( p9|) 

/ v 2d/(2-d) 

Fa d i w ( - J . (28) 
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VI. NUMERICAL RESULTS 



We carried out simulations with a simplified model: Atoms are deposited onto a one- 
dimensional surface with a deposition rate F and diffuse with a diffusion constant D until 
they meet other adatoms or a cluster of adatoms. If the encountered cluster consists of %* 
or more atoms, the atom is incorporated into this cluster, and a stable, immobile island is 
formed. Clusters of size i* or less are allowed to decay. There are no barriers for interlayer 
transport (Ehrlich-Schwoebel barriers P^| ) and no overhangs or holes are allowed (solid- 
on-solid condition). In the simulations we use a and a± as length and height unit, i.e. 
a = a± = 1. 

The squared surface width w 2 = (h 2 ) — (h) 2 as a function of time for various D/F and 
i* — 1 is shown in fig. [I] ((. . .) denotes the ensemble and spatial average). 




Ft 

FIG. 1. Squared surface width at integer and half-integer times (in units of the layer completion 
time t = 1/F) for different values of the parameter D/F = 10 5 . . . 10 10 (from top to bottom). 

For a given D/F, the oscillations in the surface width persist up to a coverage Ft, which 
increases with D/F. Beyond Ft, the crossover to kinetic roughening is observed, where w 2 
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approaches a power law t 2/3 with the cKPZ prediction of (3 = 1/3 in one dimension [fi,|l8,19 . 
Rescaling Ft by Ft = (D/F) 1 ^ 3 in fig.|2] we find an excellent collapse of the crossover regions 
for all curves of fig.[l|. This means that 5 = 1/3 within numerical accuracy, in agreement 
with pp and 7 = 1/4 |TJ. 



10' 



10 L 



10- 1 
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10' 



t/(D/F) 

FIG. 2. Curves from fig. ||, with time scaled by (D/F) 1 / 3 . 

We measured the damping time for different values of i* by determining the coverage 
Ft, where w = 0.71, 0.57 or 0.65 for i* = 1, 2, 3, respectively. It depends on I measured 
as the system size divided by the total number of nucleation events which occurred in one 
layer. Fig.^ shows that Ft ~ / 4 / 3 , independent of the values i* = 1, 2, 3, in agreement with 
the theoretical result 
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FIG. 3. Coverage Ft, at which the surface width reaches a given value (see text), as a function 
of the diffusion length I for different values i*. The straight lines are fits to the last four data 
points in each set of data. Their slopes are 1.39 ± 0.09, 1.34 ± 0.09 and 1.38 ± 0.09 for i* = 1, 2, 3, 
respectively. 

In order to check that the damping time and the layer coherence length are the appro- 
priate scales also for other quantities showing oscillations during the layer-by- layer growth 
we investigated the Bragg intensity or kinematic intensity / = {{n even — n dd) 2 ) where n even 
{ n odd) denotes the number of surface sites in even (odd) layers. It is related to the RHEED 
intensity in experiments. Fig. |] shows / at integer times, rescaled in the same way as in fig.|2]. 
Again, we find that the number of observable oscillations varies with the growth conditions 
as described by fllTf). 
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FIG. 4. Maxima of the kinematic intensity for D/F = 10 5 . . . 10 10 (right to left), and i* = 1. 
Time is scaled by (D/F) 1 / 3 . 

Finally, we carried out a finite size analysis to measure the layer coherence length I 
explicitly. As mentioned in section H, the surface does not roughen when the system size 
L is smaller than I. Then, the oscillations of the surface width w persist forever. After a 
transient time, the amplitude of the oscillations becomes stationary. We take the variance 
of the surface width w(t) during the layer completion time r = 1/F, 



as a measure of the squared amplitude of the oscillations. (. . -)[t,t+r] means the time average 
over the interval [t, t+r]. If this variance becomes equal to the ensemble fluctuations of w at 
fixed time, no oscillations can be observed. We find that A(t) approaches a stationary value 
which decreases with increasing system size. For system sizes larger than a certain value L 
A(t) is equal to the statistical fluctuations of w itself. This means that in a system of size 
L > L the oscillations can die out (or rather cannot be distinguished from noise anymore for 
long times). Therefore, L can be identified with I. According to flUf ) and ([!]) with 7 = 1/4 



A{t) 2 = (w 2 ) u+T] - (w)l t+T] , 



(29) 
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in one dimension for i* = 1 fll |, one expects L ~ {D/F) 1 / 3 . Indeed, the simulation results 
shown in fig. [5] are in excellent agreement with the theoretical prediction. 




D/F 

FIG. 5. L as a function of D/F. The fit has a slope of 0.339 ± 0.006. 
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APPENDIX 

We consider the following continuum equation : 

^ = z/V 2 /i - K{V 2 fh - AV 2 (V/i) 2 +v + Vc (30) 
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where r\ and i] c are gaussian distributed random forces with zero mean and second moment 
fa(x, t) V (x, t)) = T 5 d (x - x') 5(t - 0, (31) 



(r /c (x,t)r 7c (x ,*')) = -V V 2 <5 d (x - x ) S(t - t), (32) 



which describe the shot noise and the diffusion noise |pq| , respectively. Dimensional argu- 
ments as in section |TTT| together with (16) and (|TT|)lead to the well known expression 



V m Tt 1 « Dn (a d a ± ) 2 (33) 



for the correlator of the conserved noise |]19 |. This implies that the conserved noise dominates 



the fluctuations only on distances shorter than the typical diffusion length / [27|. As we are 
dealing with larger length scales, the conserved noise may be neglected in the following. 
The physical dimensions of the remaining parameters in (|30| ) are 



[v] = x 2 r\ [K] = xH-\ [A] = xH- l h-\ [F] = x d r l h 2 . (34) 
Comparing those of v and T one gets 

h v (t) = (F/v) d l\Ftf- d V\ (35) 



comparing those of A and T one gets [28 



h x (t) = ^ /X) d/(S+d)^ t) (A- d )/(8 +d )^ (36) 

and finally comparing those of K and T one gets 

h K (t) = (FlK) d, \Ft) {i - d)/ \ (37) 

The dimensional analysis of the linear equations, leading to ( pop and (pTf), already gives 
the right scaling behaviour of h as function of t, due to the non-renormalisation of the 
parameters u, K and T |]21j] . 

Setting h\(t\ u ) = h u (t\ u ) gives the crossover time t\ u in (|2"1|). In the same fashion, by 
setting h K (t KX ) = h x (t KX ) one gets 
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(l) =»^t) (38) 

for the crossover time from K— to A— dominated roughening. In the last equality eqns. ( [TDD 
and ([TI|) have been used. Then the crossover time agrees with the expression (P5]), consistent 
with the fact, that the i^-term and the A-term give the same result. 

Finally, one can ask for the typical times, where hx-, h\ or h v become ~ 1, i.e. the times 
which can be interpreted as the damping times, if only the corresponding term is present: 

Tt K = a ± I -y- I , Tt\ = a ± I — I , Tt v = a ± I — I . (39) 
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